*This file performs analyses using US EPA National Lakes Data
*To replicate tables and figures, please merge NSRE and NRI data by county.

*Code that uses "RegNL_Main.dta" is based on the sample that removes a random sample of "no" recreation records in order ///
*to preserve percentage of yes/no responses to water-based recreation from original NSRE data.

*Code that uses "RegNL_Placebo.dta" is based on the entire NSRE sample.

********************************************************************************************************************
*                                    Variable Definitions
********************************************************************************************************************

*Recreation 
*rec: dummy variable for water-based recreation participation
*recd: days of water-based recreation 
*lrecd: ln of days of water-based recreation + 1
*swim1,wfish1,cfish1,mboat1,visit1: participation in swimming, warmwater fishing, coldwater fishing, boating with a motor, visiting a waterbody
*tesp: outdoor team sports
*bike: bicycling
*horse: horseback riding
*picnic: picnicking
*nat: visit nature centers
*prehis: visit prehistoric sites
*hist: visit historical sites
*walk: walk for pleasure
*hike: day hiking
*wildern: visit a wilderness area
*gathermb: gather mushrooms, berries, etc.
*hunt: hunting
*siact: snow/ice activities
*mtnbike: mountain biking
*bacpac: backpacking
*devcam: developed camping
*campri: primitive camping
*bird: view or photograph birds
*fv: view or photograph fish
*wv: view/photograph other wildlife
*ov: view/photograph flowers, etc.
*scenery: view/photograph natural scenery
*driving: driving for pleasure
*mv: drive off-road
*swimp: swim in an outdoor pool

*Phosphorus and water quality
*ptl: average county-level phosphorus concentration in lakes
*streamp: average county-level phosphorus concentration in streams
*ivp0_200: upstream phosphorus concentration for 0 - 200 miles
*ivp200: upstream phosphorus concentration for 30 - 200 miles
*ivp60_200: upstream phosphorus concentration for 60 - 200 miles
*totalp1: IA Lakes phosphorus concentration
*soilloss: NRI soil loss estimate

*County info
*fips:	county fips codes
*msa: metropolitan statistical area dummy variable
*msa2:	micropolitan statistical area dummy variable
*shore:	dummy for shoreline county
*nlres: dummy for county with an EPA lake
*tmean: mean temperature
*tmeansq: mean temperature^2
*precip: mean precipitation
*zmean: mean elevation
*landcover types
*water: water
*dev: developed
*forest: forested land
*shrub: shrubland
*herb: herbaceous
*cult: planted/cultivated
*barren: barren 
*wetland: wetland

*Census Divisions
*cdiv1: northeast
*cdiv2: mid. atlantic
*cdiv3: east north central
*cdiv4: west north central
*cdiv5: south atlantic
*cdiv6: east south central
*cdiv7: west south central
*cdiv8: mountain
*cdiv9: pacific

*States:
*s2: AL
*s3: AR
*s7: CT
*s9: DE
*s10: FL
*s11: GA
*s13: IA
*s14: ID
*s15: IL
*s16: IN
*s17: KS
*s18: KY
*s19: LA
*s20: MA
*s21: MD
*s22: ME
*s23: MI
*s24: MN
*s25: MO
*s26: MS
*s27: MT
*s28: NC
*s29: ND
*s30: NE
*s31: NH
*s32: NJ
*s35: NY
*s36: OH
*s37: OK
*s38: OR
*s39: PA
*s40: RI
*s41: SC
*s42: SD
*s43: TN
*s44: TX
*s46: VA & DC
*s47: VT
*s48: WA
*s49: WI
*s50: WV

*Sociodemographics
*r1: White 
*r2: African American
*r3: American Indian or Alaska Native
*r4: Asian
*r5: Native Hawaiian or Pacific Islander
*r8: Don't Know
*r9: refused
*sex: dummy variable for male gender respondent
*age1 - age6: "age category" as labeled by NSRE
*educ1: 8th grade or less
*educ2: 9th - 11th grade
*educ3: high school graduate
*educ4: some college
*educ5: associate's degree
*educ6: bachelor's degree
*educ7: master's degree
*educ8: professional degree
*educ9: doctorate degree
*educ10: other
*educ11: don't know
*educ12: refused
*income1: <= $4,999
*income2: $5,000 to $9,999
*income3: $10,000 to $14,999
*income4: $15,000 to $19,999
*income5: $20,000 to $24,999
*income6: $25,000 to $34,999
*income7: $35,000 to $49,999
*income8: $50,000 to $74,999
*income9: $75,000 to $99,999
*income10: $100,000 to $149,999
*income11: $150,000 or more
*income12: don't know


*Define Controls
gl states "s2-s3 s7 s9-s11 s13-s32 s35-s44 s46-s49"
gl controls "msa msa2 shore nlres tmean tmeansq  precip zmean water dev forest shrub herb cult r1-r5 r8 sex age1-age6 educ1-educ9 income1-income12"

*Define file paths (corresponding to user's folders)
gl DATA "G:\Box Sync\MissingBenefits\data"
gl param "G:\Box Sync\MissingBenefits\data\ParameterEstimates"
gl FIGURES "G:\Box Sync\MissingBenefits\draft\Figures"

********************************************************************************************************************
*                                     (Table 1, Figure 1)
********************************************************************************************************************
use "$DATA\Regressions\RegNL_Main.dta", clear

*Table 1
tabstat rec recd swim1 wfish1 cfish1 mboat1 visit1 ptl, by(division)  statistics(mean, sd, n)

*Data for Figure 1
collapse ptl rec, by(fips)

********************************************************************************************************************
*                                      (Table 2)
********************************************************************************************************************

use "$DATA\Regressions\RegNL_Main.dta", clear

*Semi-Log 
regress lrecd ptl $controls $states, cluster(fips) robust 

*Semi-Log (IV)
ivreg2 lrecd $controls $states  (ptl = ivp0_200), cluster(fips) robust  

*Semi-Log (IV)
ivreg2 lrecd $controls cdiv2-cdiv9  (ptl = ivp0_200), cluster(fips) robust  

*Semi-Log (IV)
ivreg2 lrecd $controls $states  (ptl = ivp200) if ivp200>0, cluster(fips) robust  

*Semi-Log (IV)
ivreg2 lrecd $controls $states  (ptl = ivp60_200) if ivp60_200>0, cluster(fips) robust  

*Semi-Log (IV)
ivreg2 lrecd $controls $states  (ptl = ivp0_200 ivp200) if ivp200>0, cluster(fips) robust  

*Negative Binomial  
nbreg recd ptl $controls $states , cluster(fips) robust 

*GMM Count

*0-200
gmm (recd*exp(-ptl*{b1} -msa*{b3} -msa2*{b3b} -shore*{b4} -nlres*{b5} -tmean*{b6} -tmeansq*{b6b} -precip*{b7}  ///
-r1*{b8} -r2*{b8b} -r3*{b8c} -r4*{b8d} -r5*{b8e} -r8*{b8f} -sex*{b9} -age1*{b10} - age2*{b10b} -age3*{b10c} -age4*{b10d} -age5*{b10e} -age6*{b10f} ///
-educ1*{b10g} -educ2*{b10h} -educ3*{b10i} -educ4*{b10j} -educ5*{b10k} -educ6*{b10l} - educ7*{b10m} - educ8*{b10n} - educ9*{b10o}  ///
-income1*{b11b} -income2*{b11c} -income3*{b11d} -income4*{b11e} -income5*{b11f} -income6*{b11g} -income7*{b11h} -income8*{b11i} -income9*{b11j} -income10*{b11k} -income11*{b11l} -income12*{b11m} ///
-s2*{b11} -s3*{b12} -s7*{b14}  ///
-s9*{b16} -s10*{b17} -s11*{b18} -s13*{b21} -s14*{b22} -s15*{b23} -s16*{b24} -s17*{b25} -s18*{b26} -s19*{b27} -s20*{b28} -s21*{b29} -s22*{b30} -s23*{b31} -s24*{b32} ///
-s25*{b33} -s26*{b34} -s27*{b35} -s28*{b36} -s29*{b37} -s30*{b38} -s31*{b39} -s32*{b40} -s35*{b41} -s36*{b42} -s37*{b43} - s38*{b44} -s39*{b45} -s40*{b46} -s41*{b47} -s42*{b48} ///
-s43*{b49} -s44*{b50} -s46*{b51} -s47*{b52} -s48*{b53} -s49*{b54} -zmean*{b58} -water*{b59} -dev*{b60} -forest*{b61} ///
-shrub*{b62} -herb*{b63} -cult*{b64} -{const})-1),instruments(ivp0_200  ///
msa msa2 shore nlres tmean tmeansq precip  r1 r2 r3 r4 r5 r8 sex age1 age2 age3 age4 age5 age6 educ1 educ2 educ3 educ4 educ5 educ6 educ7 educ8 ///
educ9 income1 income2 income3 income4 income5 income6 income7 income8 income9 income10 income11 income12 ///
s2 s3  s7 s9 s10 s11 s13 s14 s15 s16 s17 s18 s19 s20 s21 s22 s23 s24 s25 ///
s26 s27 s28 s29 s30 s31 s32 s35 s36 s37 s38 s39 s40 s41 s42 s43 s44 s46 s47 s48 s49 zmean water dev forest shrub herb cult ) wmatrix(cluster fips) twostep nolog 

*0-200, census division FE
gmm (recd*exp(-ptl*{b1} -msa*{b3} -msa2*{b3b} -shore*{b4} -nlres*{b5} -tmean*{b6} -tmeansq*{b6b} -precip*{b7}  ///
-r1*{b8} -r2*{b8b} -r3*{b8c} -r4*{b8d} -r5*{b8e} -r8*{b8f} -sex*{b9} -age1*{b10} - age2*{b10b} -age3*{b10c} -age4*{b10d} -age5*{b10e} -age6*{b10f} ///
-educ1*{b10g} -educ2*{b10h} -educ3*{b10i} -educ4*{b10j} -educ5*{b10k} -educ6*{b10l} - educ7*{b10m} - educ8*{b10n} - educ9*{b10o}  ///
-income1*{b11b} -income2*{b11c} -income3*{b11d} -income4*{b11e} -income5*{b11f} -income6*{b11g} -income7*{b11h} -income8*{b11i} -income9*{b11j} -income10*{b11k} -income11*{b11l} -income12*{b11m} ///
-cdiv2*{b11} -cdiv3*{b12} -cdiv4*{b14}  ///
-cdiv5*{b16} -cdiv6*{b17} -cdiv7*{b18} -cdiv8*{b21} -cdiv9*{b22} -zmean*{b58} -water*{b59} -dev*{b60} -forest*{b61} ///
-shrub*{b62} -herb*{b63} -cult*{b64} -{const})-1),instruments(ivp0_200  ///
msa msa2 shore nlres tmean tmeansq precip r1 r2 r3 r4 r5 r8 sex age1 age2 age3 age4 age5 age6 educ1 educ2 educ3 educ4 educ5 educ6 educ7 educ8 ///
educ9 income1 income2 income3 income4 income5 income6 income7 income8 income9 income10 income11 income12 ///
cdiv2 cdiv3 cdiv4 cdiv5 cdiv6 cdiv7 cdiv8 cdiv9 ///
zmean water dev forest shrub herb cult ) wmatrix(cluster fips) twostep nolog 

*30-200
gmm (recd*exp(-ptl*{b1} -msa*{b3} -msa2*{b3b} -shore*{b4} -nlres*{b5} -tmean*{b6} -tmeansq*{b6b} -precip*{b7}  ///
-r1*{b8} -sex*{b9} -age1*{b10} - age2*{b10b} -age3*{b10c} -age4*{b10d} -age5*{b10e} -age6*{b10f} ///
-educ1*{b10g} -educ2*{b10h} -educ3*{b10i} -educ4*{b10j} -educ5*{b10k} -educ6*{b10l} - educ7*{b10m} - educ8*{b10n} - educ9*{b10o}  ///
-income1*{b11b} -income2*{b11c} -income3*{b11d} -income4*{b11e} -income5*{b11f} -income6*{b11g} -income7*{b11h} -income8*{b11i} -income9*{b11j} -income10*{b11k} -income11*{b11l} -income12*{b11m} ///
-s2*{b11} -s3*{b12} -s7*{b14}  ///
-s9*{b16} -s10*{b17} -s11*{b18} -s13*{b21} -s14*{b22} -s15*{b23} -s16*{b24} -s17*{b25} -s18*{b26} -s19*{b27} -s20*{b28} -s21*{b29} -s22*{b30} -s23*{b31} -s24*{b32} ///
-s25*{b33} -s26*{b34} -s27*{b35} -s28*{b36} -s29*{b37} -s30*{b38} -s31*{b39} -s32*{b40} -s35*{b41} -s36*{b42} -s37*{b43} - s38*{b44} -s39*{b45} -s40*{b46} -s41*{b47} -s42*{b48} ///
-s43*{b49} -s44*{b50} -s46*{b51} -s47*{b52} -s48*{b53} -s49*{b54} -zmean*{b58} -water*{b59} -dev*{b60} -forest*{b61} ///
-shrub*{b62} -herb*{b63} -cult*{b64} -{const})-1) if ivp200>0,instruments(ivp200   ///
msa msa2 shore nlres tmean tmeansq precip r1 sex age1 age2 age3 age4 age5 age6 educ1 educ2 educ3 educ4 educ5 educ6 educ7 educ8 ///
educ9 income1 income2 income3 income4 income5 income6 income7 income8 income9 income10 income11 income12 ///
s2 s3  s7 s9 s10 s11 s13 s14 s15 s16 s17 s18 s19 s20 s21 s22 s23 s24 s25 ///
s26 s27 s28 s29 s30 s31 s32 s35 s36 s37 s38 s39 s40 s41 s42 s43 s44 s46 s47 s48 s49 zmean water dev forest shrub herb cult ) wmatrix(cluster fips) twostep nolog 

*60-200
gmm (recd*exp(-ptl*{b1} -msa*{b3} -msa2*{b3b} -shore*{b4} -nlres*{b5} -tmean*{b6} -tmeansq*{b6b} -precip*{b7}  ///
-r1*{b8} -sex*{b9} -age1*{b10} - age2*{b10b} -age3*{b10c} -age4*{b10d} -age5*{b10e} -age6*{b10f} ///
-educ1*{b10g} -educ2*{b10h} -educ3*{b10i} -educ4*{b10j} -educ5*{b10k} -educ6*{b10l} - educ7*{b10m} - educ8*{b10n} - educ9*{b10o}  ///
-income1*{b11b} -income2*{b11c} -income3*{b11d} -income4*{b11e} -income5*{b11f} -income6*{b11g} -income7*{b11h} -income8*{b11i} -income9*{b11j} -income10*{b11k} -income11*{b11l} -income12*{b11m} ///
-s2*{b11} -s3*{b12} -s7*{b14}  ///
-s9*{b16} -s10*{b17} -s11*{b18} -s13*{b21} -s14*{b22} -s15*{b23} -s16*{b24} -s17*{b25} -s18*{b26} -s19*{b27} -s20*{b28} -s21*{b29} -s22*{b30} -s23*{b31} -s24*{b32} ///
-s25*{b33} -s26*{b34} -s27*{b35} -s28*{b36} -s29*{b37} -s30*{b38} -s31*{b39} -s32*{b40} -s35*{b41} -s36*{b42} -s37*{b43} - s38*{b44} -s39*{b45}  -s41*{b47} -s42*{b48} ///
-s43*{b49} -s44*{b50} -s46*{b51} -s47*{b52} -s48*{b53} -s49*{b54} -zmean*{b58} -water*{b59} -dev*{b60} -forest*{b61} ///
-shrub*{b62} -herb*{b63} -cult*{b64} -{const})-1) if ivp60_200>0,instruments(ivp60_200   ///
msa msa2 shore nlres tmean tmeansq precip r1 sex age1 age2 age3 age4 age5 age6 educ1 educ2 educ3 educ4 educ5 educ6 educ7 educ8 ///
educ9 income1 income2 income3 income4 income5 income6 income7 income8 income9 income10 income11 income12 ///
s2 s3  s7 s9 s10 s11 s13 s14 s15 s16 s17 s18 s19 s20 s21 s22 s23 s24 s25 ///
s26 s27 s28 s29 s30 s31 s32 s35 s36 s37 s38 s39 s41 s42 s43 s44 s46 s47 s48 s49 zmean water dev forest shrub herb cult ) wmatrix(cluster fips) twostep nolog 

*all IV
gmm (recd*exp(-ptl*{b1} -msa*{b3} -msa2*{b3b} -shore*{b4} -nlres*{b5} -tmean*{b6} -tmeansq*{b6b} -precip*{b7}  ///
-r1*{b8} -sex*{b9} -age1*{b10} - age2*{b10b} -age3*{b10c} -age4*{b10d} -age5*{b10e} -age6*{b10f} ///
-educ1*{b10g} -educ2*{b10h} -educ3*{b10i} -educ4*{b10j} -educ5*{b10k} -educ6*{b10l} - educ7*{b10m} - educ8*{b10n} - educ9*{b10o}  ///
-income1*{b11b} -income2*{b11c} -income3*{b11d} -income4*{b11e} -income5*{b11f} -income6*{b11g} -income7*{b11h} -income8*{b11i} -income9*{b11j} -income10*{b11k} -income11*{b11l} -income12*{b11m} ///
-s2*{b11} -s3*{b12} -s7*{b14}  ///
-s9*{b16} -s10*{b17} -s11*{b18} -s13*{b21} -s14*{b22} -s15*{b23} -s16*{b24} -s17*{b25} -s18*{b26} -s19*{b27} -s20*{b28} -s21*{b29} -s22*{b30} -s23*{b31} -s24*{b32} ///
-s25*{b33} -s26*{b34} -s27*{b35} -s28*{b36} -s29*{b37} -s30*{b38} -s31*{b39} -s32*{b40} -s35*{b41} -s36*{b42} -s37*{b43} - s38*{b44} -s39*{b45} -s40*{b46} -s41*{b47} -s42*{b48} ///
-s43*{b49} -s44*{b50} -s46*{b51} -s47*{b52} -s48*{b53} -s49*{b54} -zmean*{b58} -water*{b59} -dev*{b60} -forest*{b61} ///
-shrub*{b62} -herb*{b63} -cult*{b64} -{const})-1) if ivp200>0 & ivp0_200>0,instruments(ivp0_200 ivp200   ///
msa msa2 shore nlres tmean tmeansq precip  r1 sex age1 age2 age3 age4 age5 age6 educ1 educ2 educ3 educ4 educ5 educ6 educ7 educ8 ///
educ9 income1 income2 income3 income4 income5 income6 income7 income8 income9 income10 income11 income12 ///
s2 s3  s7 s9 s10 s11 s13 s14 s15 s16 s17 s18 s19 s20 s21 s22 s23 s24 s25 ///
s26 s27 s28 s29 s30 s31 s32 s35 s36 s37 s38 s39 s40 s41 s42 s43 s44 s46 s47 s48 s49 zmean water dev forest shrub herb cult ) wmatrix(cluster fips) twostep nolog 
estat overid


********************************************************************************************************************
*                                      (Table 4 - Reduced Form and First Stage)
********************************************************************************************************************

use "$DATA\Regressions\RegNL_Main.dta", clear

*Reduced-Form
*Semi-Log 
regress lrecd ivp0_200 $controls $states , cluster(fips) robust 

*Semi-Log 
regress lrecd ivp200  $controls $states if ivp200>0, cluster(fips) robust 
 
*Semi-Log 
regress lrecd ivp60_200 $controls $states if ivp60_200>0 , cluster(fips) robust 
 
*First Stage
*Semi-Log 
regress ptl ivp0_200 $controls $states , cluster(fips) robust 
 
*Semi-Log 
regress ptl ivp200  $controls $states if ivp200>0, cluster(fips) robust 
 
*Semi-Log 
regress ptl ivp60_200 $controls $states if ivp60_200>0 , cluster(fips) robust 


********************************************************************************************************************
*                                      (Tables 5,6,7; Figures 2,3 - Measurement Error)
********************************************************************************************************************
********************************************************************
* Table 5 Coefficient estimates for models with single regressors  *
********************************************************************
use "$DATA\Regressions\RegNL_Main.dta", clear

*Semi-Log 
regress lrecd ptl , cluster(fips) robust 

*Semi-Log (IV)
ivreg2 lrecd  (ptl  = ivp0_200), cluster(fips) robust 

*Negative Binomial  
nbreg recd ptl , cluster(fips) robust 

*GMM Count
*Clustered at county level 
gmm (recd*exp(-ptl*{b1} -{const})-1),instruments(ivp0_200) wmatrix(cluster fips) twostep nolog 
estat overid

********************************************************************
* Table 5 - Calculating Attenuation Factors						   *
********************************************************************

use "$DATA\Regressions\RegNL_Placebo.dta", clear
collapse ptl totalp1 ivp0_200 $controls , by(fips)

keep if !missing(totalp1) & !missing(ptl)
gen perr1 = ptl-totalp1	

*Export to Excel to Calculate Variance and Covariance


********************************************************************
* Table 6,7 - M.E. vs. Local Concentrations  and IV				   *
********************************************************************

*No Controls
regress perr1 totalp1   , robust

regress perr1 ptl  , robust

regress perr1 ivp0_200  , robust

*Controls
regress perr1 totalp1  $controls , robust

regress perr1 ptl  $controls , robust

regress perr1 ivp0_200  $controls , robust


********************************************************************
* Figures 2,3													   *
********************************************************************

*Figure 2
histogram perr1, bin(20) frequency fcolor(ltbluishgray) lcolor(black) lpattern(solid) graphregion(fcolor(white)) ///
kdensity kdenopts(lcolor(red) lpattern(dash)) xtitle(Measurement Error (mg/l)) 

*Figure 3
*Measurement Error vs. National Lakes
twoway (scatter perr1 ptl, mcolor(navy)) (lfit perr1 ptl, lcolor(red) fcolor(white)), ytitle(Measurement Error (mg/l)) ///
xtitle(US EPA National Lakes Concentrations (mg/l)) graphregion(fcolor(white) lpattern(blank)) yscale(range(-.4 .6)) ylabel(-.4(.2).6) ///
legend(off)

*Measurement Error vs. True Concentration (IA Lakes)
twoway (scatter perr1 totalp1, mcolor(navy)) (lfit perr1 totalp1, lcolor(red) fcolor(white)), ytitle(Measurement Error (mg/l)) ///
xtitle(Iowa Lakes Concentrations (mg/l)) graphregion(fcolor(white) lpattern(blank)) yscale(range(-.4 .6)) ylabel(-.4(.2).6) ///
legend(off)

*Measurement Error vs. Upstream IV
twoway (scatter perr1 ivp0_200, mcolor(navy)) (lfit perr1 ivp0_200, lcolor(red) fcolor(white)), ytitle(Measurement Error (mg/l)) ///
xtitle(Upstream IV Concentrations (mg/l)) graphregion(fcolor(white) lpattern(blank)) yscale(range(-.4 .6)) ylabel(-.4(.2).6) ///
legend(off)


********************************************************************************************************************
*                                      (Appendix Tables 1,3 - Placebo Tests, with and without weights)
********************************************************************************************************************
use "$DATA\Regressions\RegNL_Main.dta", clear


foreach act in "tesp" "bike" "horse" "picnic" "nat" "prehis" "hist" "walk" "hike" "wildern" "gathermb" "hunt" "siact" "mtnbike" "bacpac" "devcam" "campri" "bird" "fv" "wv" "ov" ///
"scenery" "driving" "mv" "swimp"    {
gen `act'yn= 0 if `act'==2
replace `act'yn=1 if `act'==1
}
*
*Columns 1,2
*Probit
ivprobit rec $controls $states  (ptl  = ivp0_200 ), cluster(fips) robust 

ivprobit rec $controls tespyn-swimpyn $states  (ptl  = ivp0_200 ), vce(cluster fips)

*Columns 1,2
ivprobit rec $controls $states  (ptl  =  ivp0_200 ) [pw=weight], cluster(fips) robust 

*Probit
ivprobit rec $controls tespyn-swimpyn $states  (ptl  =  ivp0_200 ) [pw=weight], vce(cluster fips) 



use "$DATA\Regressions\RegNL_Placebo.dta", clear

foreach act in "tesp" "bike" "horse" "picnic" "nat" "prehis" "hist" "walk" "hike" "wildern" "gathermb" "hunt" "siact" "mtnbike" "bacpac" "devcam" "campri" "bird" "fv" "wv" "ov" ///
"scenery" "driving" "mv" "swimp"    {
gen `act'yn= 0 if `act'==2
replace `act'yn=1 if `act'==1
}
*

*Columns 3, 4
preserve
*Keep records if individual answered Y/N to all swimming, fishing, boating, visitation questions
keep if (swim==1 | swim==2) & (cfish==1 | cfish==2) & (wfish==1 | wfish==2) & (mboat==1 | mboat==2) & (visit==1 | visit==2)

*Create field for recreating or not
gen rec=0
replace rec=1 if swim==1 | cfish==1 | wfish==1 | mboat==1 | visit==1

*Probit
ivprobit rec $controls $states  (ptl  = ivp0_200 ) , robust cluster(fips)

*Probit
ivprobit rec $controls $states (ptl  = ivp0_200 ) [pw=weight], robust cluster(fips)

*Probit
ivprobit rec $controls $states tespyn-swimpyn  (ptl  = ivp0_200 ) , vce(cluster fips)

*Probit
ivprobit rec $controls $states tespyn-swimpyn (ptl  = ivp0_200 ) [pw=weight], vce(cluster fips)

restore

gl ctrl_bike "tespyn horseyn-swimpyn"
gl ctrl_hike "tespyn-walkyn wildernyn-swimpyn"
gl ctrl_siact "tespyn-huntyn mtnbikeyn-swimpyn"
gl ctrl_swimp "tespyn-mvyn"
gl ctrl_tesp "bikeyn-swimpyn"
gl ctrl_horse "tespyn-bikeyn picnicyn-swimpyn"
gl ctrl_picnic "tespyn-horseyn natyn-swimpyn"
gl ctrl_nat "tespyn-picnicyn prehisyn-swimpyn"
gl ctrl_prehis "tespyn-natyn histyn-swimpyn"
gl ctrl_hist "tespyn-prehisyn walkyn-swimpyn"
gl ctrl_walk "tespyn-histyn wildernyn-swimpyn"
gl ctrl_wildern "tespyn-walkyn gathermbyn-swimpyn"
gl ctrl_gathermb "tespyn-wildernyn huntyn-swimpyn"
gl ctrl_hunt "tespyn-gathermbyn mtnbikeyn-swimpyn"
gl ctrl_mtnbike "tespyn-huntyn bacpacyn-swimpyn"
gl ctrl_bacpac "tespyn-mtnbikeyn devcamyn-swimpyn"
gl ctrl_devcam "tespyn-bacpacyn campriyn-swimpyn"
gl ctrl_campri "tespyn-devcamyn birdyn-swimpyn"
gl ctrl_bird "tespyn-campriyn fvyn-swimpyn"
gl ctrl_fv "tespyn-birdyn wvyn-swimpyn"
gl ctrl_wv "tespyn-fvyn ovyn-swimpyn"
gl ctrl_ov "tespyn-wvyn sceneryyn-swimpyn"
gl ctrl_scenery "tespyn-ovyn drivingyn-swimpyn"
gl ctrl_driving "tespyn-sceneryyn mvyn-swimpyn"
gl ctrl_mv "tespyn-drivingyn swimpyn"

*NL
foreach act in "bike" "hike" "siact"  "swimp" {
ivprobit `act'yn $controls $states (ptl  =  ivp0_200 ) ,vce(cluster fips)
}
*

foreach act in "swim1" "cfish1" "wfish1" "mboat1" "visit1" {
ivprobit `act' $controls $states  (ptl  = ivp0_200 ) ,vce(cluster fips)
}
*
foreach act in "tesp" "horse" "picnic" "nat" "prehis" "hist" "walk" "wildern" "gathermb" "hunt" "mtnbike" "bacpac" "devcam" "campri" "bird" "fv" "wv" "ov" ///
"scenery" "driving" "mv" {
ivprobit `act'yn $controls $states (ptl = ivp0_200  ),vce(cluster fips)
}
*

*Bike
ivprobit bikeyn $controls $states $ctrl_bike (ptl  =  ivp0_200 ) ,vce(cluster fips)

*Hike
ivprobit hikeyn $controls $states $ctrl_hike (ptl  =  ivp0_200 ) ,vce(cluster fips)

*Siact
ivprobit siactyn $controls $states $ctrl_siact (ptl  =  ivp0_200 ) ,vce(cluster fips)

*Swimp
ivprobit swimpyn $controls $states $ctrl_swimp (ptl  =  ivp0_200 ) ,vce(cluster fips)


foreach act in "swim1" "cfish1" "wfish1" "mboat1" "visit1" {
ivprobit `act' $controls $states tespyn-swimpyn  (ptl  = ivp0_200 ) ,vce(cluster fips)
}
*

*tesp
ivprobit tespyn $controls $states $ctrl_tesp (ptl  =  ivp0_200 ) ,vce(cluster fips)

*horse
ivprobit horseyn $controls $states $ctrl_horse (ptl  =  ivp0_200 ) ,vce(cluster fips)

*picnic
ivprobit picnicyn $controls $states $ctrl_picnic (ptl  =  ivp0_200 ) ,vce(cluster fips)

*nat
ivprobit natyn $controls $states $ctrl_nat (ptl  =  ivp0_200 ) ,vce(cluster fips)

*prehis
ivprobit prehisyn $controls $states $ctrl_prehis (ptl  =  ivp0_200 ) ,vce(cluster fips)

*hist
ivprobit histyn $controls $states $ctrl_hist (ptl  =  ivp0_200 ) ,vce(cluster fips)

*walk
ivprobit walkyn $controls $states $ctrl_walk (ptl  =  ivp0_200 ) ,vce(cluster fips)

*wildern
ivprobit wildernyn $controls $states $ctrl_wildern (ptl  =  ivp0_200 ) ,vce(cluster fips)

*gathermb
ivprobit gathermbyn $controls $states $ctrl_gathermb (ptl  =  ivp0_200 ) ,vce(cluster fips)

*hunt
ivprobit huntyn $controls $states $ctrl_hunt (ptl  =  ivp0_200 ) ,vce(cluster fips)

*mtnbike
ivprobit mtnbikeyn $controls $states $ctrl_mtnbike (ptl  =  ivp0_200 ) ,vce(cluster fips)

*bacpac
ivprobit bacpacyn $controls $states $ctrl_bacpac (ptl  =  ivp0_200 ) ,vce(cluster fips)

*devcam
ivprobit devcamyn $controls $states $ctrl_devcam (ptl  =  ivp0_200 ) ,vce(cluster fips)

*campri
ivprobit campriyn $controls $states $ctrl_campri (ptl  =  ivp0_200 ) ,vce(cluster fips)

*bird
ivprobit birdyn $controls $states $ctrl_bird (ptl  =  ivp0_200 ) ,vce(cluster fips)

*fv
ivprobit fvyn $controls $states $ctrl_fv (ptl  =  ivp0_200 ) ,vce(cluster fips)

*wv
ivprobit wvyn $controls $states $ctrl_wv (ptl  =  ivp0_200 ) ,vce(cluster fips)

*ov
ivprobit ovyn $controls $states $ctrl_ov (ptl  =  ivp0_200 ) ,vce(cluster fips)

*scenery
ivprobit sceneryyn $controls $states $ctrl_scenery (ptl  =  ivp0_200 ) ,vce(cluster fips)

*driving
ivprobit drivingyn $controls $states $ctrl_driving (ptl  =  ivp0_200 ) ,vce(cluster fips)

*mv
ivprobit mvyn $controls $states $ctrl_mv (ptl  =  ivp0_200 ) ,vce(cluster fips)


*weights
foreach act in "bike" "hike" "siact"  "swimp" {
ivprobit `act'yn $controls $states (ptl  = ivp0_200 ) [pw=weight] ,vce(cluster fips) 
}
*

foreach act in "swim1" "cfish1" "wfish1" "mboat1" "visit1" {
ivprobit `act' $controls $states (ptl  = ivp0_200 ) [pw=weight] ,vce(cluster fips)
}
*
foreach act in "tesp" "horse" "picnic" "nat" "prehis" "hist" "walk" "wildern" "gathermb" "hunt" "mtnbike" "bacpac" "devcam" "campri" "bird" "fv" "wv" "ov" ///
"scenery" "driving" "mv" {
ivprobit `act'yn $controls $states (ptl = ivp0_200  ) [pw=weight] ,vce(cluster fips)
}
*

foreach act in "swim1" "cfish1" "wfish1" "mboat1" "visit1" {
ivprobit `act' $controls $states tespyn-swimpyn  (ptl  = ivp0_200 ) [pw=weight] ,vce(cluster fips)
}
*
*Bike
ivprobit bikeyn $controls $states $ctrl_bike  (ptl  =  ivp0_200 ) [pw=weight] ,vce(cluster fips)

*Hike
ivprobit hikeyn $controls $states $ctrl_hike  (ptl  =  ivp0_200 ) [pw=weight] ,vce(cluster fips)

*Siact
ivprobit siactyn $controls $states $ctrl_siact  (ptl  =  ivp0_200 ) [pw=weight] ,vce(cluster fips)

*Swimp
ivprobit swimpyn $controls $states $ctrl_swimp  (ptl  =  ivp0_200 ) [pw=weight] ,vce(cluster fips)

*tesp
ivprobit tespyn $controls $states $ctrl_tesp  (ptl  =  ivp0_200 ) [pw=weight] ,vce(cluster fips)

*horse
ivprobit horseyn $controls $states $ctrl_horse  (ptl  =  ivp0_200 ) [pw=weight] ,vce(cluster fips)

*picnic
ivprobit picnicyn $controls $states $ctrl_picnic  (ptl  =  ivp0_200 ) [pw=weight] ,vce(cluster fips)

*nat
ivprobit natyn $controls $states $ctrl_nat  (ptl  =  ivp0_200 ) [pw=weight] ,vce(cluster fips)

*prehis
ivprobit prehisyn $controls $states $ctrl_prehis  (ptl  =  ivp0_200 ) [pw=weight] ,vce(cluster fips)

*hist
ivprobit histyn $controls $states $ctrl_hist (ptl  =  ivp0_200 ) [pw=weight] ,vce(cluster fips)

*walk
ivprobit walkyn $controls $states $ctrl_walk  (ptl  =  ivp0_200 ) [pw=weight] ,vce(cluster fips)

*wildern
ivprobit wildernyn $controls $states $ctrl_wildern  (ptl  =  ivp0_200 ) [pw=weight] ,vce(cluster fips)

*gathermb
ivprobit gathermbyn $controls $states $ctrl_gathermb  (ptl  =  ivp0_200 ) [pw=weight] ,vce(cluster fips)

*hunt
ivprobit huntyn $controls $states $ctrl_hunt (ptl  =  ivp0_200 ) [pw=weight] ,vce(cluster fips)

*mtnbike
ivprobit mtnbikeyn $controls $states $ctrl_mtnbike  (ptl  =  ivp0_200 ) [pw=weight] ,vce(cluster fips)

*bacpac
ivprobit bacpacyn $controls $states $ctrl_bacpac  (ptl  =  ivp0_200 ) [pw=weight] ,vce(cluster fips)

*devcam
ivprobit devcamyn $controls $states $ctrl_devcam  (ptl  =  ivp0_200 ) [pw=weight] ,vce(cluster fips)

*campri
ivprobit campriyn $controls $states $ctrl_campri  (ptl  =  ivp0_200 ) [pw=weight] ,vce(cluster fips)

*bird
ivprobit birdyn $controls $states $ctrl_bird  (ptl  =  ivp0_200 ) [pw=weight] ,vce(cluster fips)

*fv
ivprobit fvyn $controls $states $ctrl_fv  (ptl  =  ivp0_200 ) [pw=weight] ,vce(cluster fips)

*wv
ivprobit wvyn $controls $states $ctrl_wv (ptl  =  ivp0_200 ) [pw=weight] ,vce(cluster fips)

*ov
ivprobit ovyn $controls $states $ctrl_ov  (ptl  =  ivp0_200 ) [pw=weight] ,vce(cluster fips)

*scenery
ivprobit sceneryyn $controls $states $ctrl_scenery  (ptl  =  ivp0_200 ) [pw=weight] ,vce(cluster fips)

*driving
ivprobit drivingyn $controls $states $ctrl_driving  (ptl  =  ivp0_200 ) [pw=weight],vce(cluster fips)

*mv
ivprobit mvyn $controls $states $ctrl_mv  (ptl  =  ivp0_200 ) [pw=weight] ,vce(cluster fips)


********************************************************************************************************************
*                                      (Appendix Table 5 - Soil Erosion)
********************************************************************************************************************
use "$DATA\Regressions\RegNL_Main.dta", clear
keep if recd>=1
gl controlsoil "income1-income12 age1-age6 sex educ1-educ9 forest"

*Semi-Log 
regress lrecd soilloss, cluster(fips) robust 

*Semi-Log 
regress lrecd soilloss $controlsoil, cluster(fips) robust 

*Semi-Log
regress lrecd soilloss $controls $states, cluster(fips) robust 
  
*Semi-Log (IV)
ivreg2 lrecd (soilloss  = ivp0_200 ), cluster(fips) robust  

*Semi-Log (IV)
ivreg2 lrecd $controlsoil (soilloss  = ivp0_200 ), cluster(fips) robust partial($controlsoil)
 
*Semi-Log (IV)
ivreg2 lrecd $controls $states (soilloss  = ivp0_200 ), cluster(fips) robust partial($controls $states) 

*Negative Binomial  
nbreg recd soilloss  , cluster(fips) robust 

*Negative Binomial  
nbreg recd soilloss $controlsoil , cluster(fips) robust 

*Negative Binomial  
nbreg recd soilloss $controls $states , cluster(fips) robust 

*GMM Count
gmm (recd*exp(-soilloss*{b1} -{const})-1),instruments(ivp0_200) wmatrix(cluster fips) twostep nolog 

*GMM Count
gmm (recd*exp(-soilloss*{b1} -sex*{b9} -age1*{b10} - age2*{b10b} -age3*{b10c} -age4*{b10d} -age5*{b10e} -age6*{b10f} ///
-educ1*{b10g} -educ2*{b10h} -educ3*{b10i} -educ4*{b10j} -educ5*{b10k} -educ6*{b10l} - educ7*{b10m} - educ8*{b10n} - educ9*{b10o}  ///
-income1*{b11b} -income2*{b11c} -income3*{b11d} -income4*{b11e} -income5*{b11f} -income6*{b11g} -income7*{b11h} -income8*{b11i} -income9*{b11j} -income10*{b11k} -income11*{b11l} -income12*{b11m} ///
-forest*{b61} ///
-{const})-1),instruments(ivp0_200  ///
sex age1 age2 age3 age4 age5 age6 educ1 educ2 educ3 educ4 educ5 educ6 educ7 educ8 ///
educ9  income1 income2 income3 income4 income5 income6 income7 income8 income9 income10 income11 income12 ///
forest  ) wmatrix(cluster fips) twostep nolog 

*GMM Count
gmm (recd*exp(-soilloss*{b1} -msa*{b3} -msa2*{b3b} -shore*{b4} -nlres*{b5} -tmean*{b6} -tmeansq*{b6b} -precip*{b7}  ///
-r1*{b8} -r2*{b8b} -r3*{b8c} -r4*{b8d} -r5*{b8e} -r8*{b8f} -sex*{b9} -age1*{b10} - age2*{b10b} -age3*{b10c} -age4*{b10d} -age5*{b10e} -age6*{b10f} ///
-educ1*{b10g} -educ2*{b10h} -educ3*{b10i} -educ4*{b10j} -educ5*{b10k} -educ6*{b10l} - educ7*{b10m} - educ8*{b10n} - educ9*{b10o}  ///
-income1*{b11b} -income2*{b11c} -income3*{b11d} -income4*{b11e} -income5*{b11f} -income6*{b11g} -income7*{b11h} -income8*{b11i} -income9*{b11j} -income10*{b11k} -income11*{b11l} -income12*{b11m} ///
-s2*{b11} -s3*{b12} -s7*{b14}  ///
-s9*{b16} -s10*{b17} -s11*{b18} -s13*{b21} -s14*{b22} -s15*{b23} -s16*{b24} -s17*{b25} -s18*{b26} -s19*{b27} -s20*{b28} -s21*{b29} -s22*{b30} -s23*{b31} -s24*{b32} ///
-s25*{b33} -s26*{b34} -s27*{b35} -s28*{b36} -s29*{b37} -s30*{b38} -s31*{b39} -s32*{b40} -s35*{b41} -s36*{b42} -s37*{b43} - s38*{b44} -s39*{b45} -s40*{b46} -s41*{b47} -s42*{b48} ///
-s43*{b49} -s44*{b50} -s46*{b51} -s47*{b52} -s48*{b53} -s49*{b54} -zmean*{b58} -water*{b59} -dev*{b60} -forest*{b61} ///
-shrub*{b62} -herb*{b63} -cult*{b64} -{const})-1),instruments(ivp0_200  ///
msa msa2 shore nlres tmean tmeansq precip  r1 r2 r3 r4 r5 r8 sex age1 age2 age3 age4 age5 age6 educ1 educ2 educ3 educ4 educ5 educ6 educ7 educ8 ///
educ9  income1 income2 income3 income4 income5 income6 income7 income8 income9 income10 income11 income12 ///
s2 s3  s7 s9 s10 s11 s13 s14 s15 s16 s17 s18 s19 s20 s21 s22 s23 s24 s25 ///
s26 s27 s28 s29 s30 s31 s32 s35 s36 s37 s38 s39 s40 s41 s42 s43 s44 s46 s47 s48 s49 zmean water dev forest shrub herb cult ) wmatrix(cluster fips) twostep nolog 


********************************************************************************************************************
*                                     Other Findings Discussed
********************************************************************************************************************

*Streams

use "$DATA\Regressions\RegNL_Main.dta", clear
*Semi-Log 
regress lrecd streamp $controls $states, cluster(fips) robust 

*Semi-Log (IV)
ivreg2 lrecd $controls $states  (streamp = ivp0_200), cluster(fips) robust  

*Semi-Log (IV)
ivreg2 lrecd $controls cdiv2-cdiv9  (streamp = ivp0_200), cluster(fips) robust  

*Semi-Log (IV)
ivreg2 lrecd $controls $states  (streamp = ivp200) if ivp200>0, cluster(fips) robust  

*Semi-Log (IV)
ivreg2 lrecd $controls $states  (streamp = ivp60_200) if ivp60_200>0, cluster(fips) robust  

*Semi-Log (IV)
ivreg2 lrecd $controls $states  (streamp = ivp0_200 ivp200) if ivp200>0, cluster(fips) robust  


*Average Marginal Effects
use "$DATA\Regressions\RegNL_Main.dta", clear

*Negative Binomial  
nbreg recd ptl $controls $states , cluster(fips) robust
margins, dydx(ptl)

*GMM Count
gmm (recd*exp(-ptl*{b1} -msa*{b3} -msa2*{b3b} -shore*{b4} -nlres*{b5} -tmean*{b6} -tmeansq*{b6b} -precip*{b7}  ///
-r1*{b8} -r2*{b8b} -r3*{b8c} -r4*{b8d} -r5*{b8e} -r8*{b8f} -sex*{b9} -age1*{b10} - age2*{b10b} -age3*{b10c} -age4*{b10d} -age5*{b10e} -age6*{b10f} ///
-educ1*{b10g} -educ2*{b10h} -educ3*{b10i} -educ4*{b10j} -educ5*{b10k} -educ6*{b10l} - educ7*{b10m} - educ8*{b10n} - educ9*{b10o}  ///
-income1*{b11b} -income2*{b11c} -income3*{b11d} -income4*{b11e} -income5*{b11f} -income6*{b11g} -income7*{b11h} -income8*{b11i} -income9*{b11j} -income10*{b11k} -income11*{b11l} -income12*{b11m} ///
-s2*{b11} -s3*{b12} -s7*{b14}  ///
-s9*{b16} -s10*{b17} -s11*{b18} -s13*{b21} -s14*{b22} -s15*{b23} -s16*{b24} -s17*{b25} -s18*{b26} -s19*{b27} -s20*{b28} -s21*{b29} -s22*{b30} -s23*{b31} -s24*{b32} ///
-s25*{b33} -s26*{b34} -s27*{b35} -s28*{b36} -s29*{b37} -s30*{b38} -s31*{b39} -s32*{b40} -s35*{b41} -s36*{b42} -s37*{b43} - s38*{b44} -s39*{b45} -s40*{b46} -s41*{b47} -s42*{b48} ///
-s43*{b49} -s44*{b50} -s46*{b51} -s47*{b52} -s48*{b53} -s49*{b54} -zmean*{b58} -water*{b59} -dev*{b60} -forest*{b61} ///
-shrub*{b62} -herb*{b63} -cult*{b64} -{const})-1),instruments(ivp0_200  ///
msa msa2 shore nlres tmean tmeansq precip  r1 r2 r3 r4 r5 r8 sex age1 age2 age3 age4 age5 age6 educ1 educ2 educ3 educ4 educ5 educ6 educ7 educ8 ///
educ9 income1 income2 income3 income4 income5 income6 income7 income8 income9 income10 income11 income12 ///
s2 s3  s7 s9 s10 s11 s13 s14 s15 s16 s17 s18 s19 s20 s21 s22 s23 s24 s25 ///
s26 s27 s28 s29 s30 s31 s32 s35 s36 s37 s38 s39 s40 s41 s42 s43 s44 s46 s47 s48 s49 zmean water dev forest shrub herb cult ) wmatrix(cluster fips) twostep nolog 

matrix est = e(b)

forvalues i = 1/90{
gen est`i' = est[1,`i']
}
*
local est1 = est[1,1]

gen ptlb = est1*ptl
gen msab = est2* msa
gen msa2b = est3* msa2
gen shoreb = est4* shore
gen nlresb = est5* nlres
gen tmeanb = est6* tmean
gen tmeansqb = est7* tmeansq
gen precipb = est8* precip
gen r1b = est9* r1
gen r2b = est10* r2
gen r3b = est11* r3
gen r4b = est12* r4
gen r5b = est13* r5
gen r8b = est14* r8
gen sexb = est15* sex
gen age1b = est16* age1
gen age2b = est17* age2
gen age3b = est18* age3
gen age4b = est19* age4
gen age5b = est20* age5
gen age6b = est21* age6
gen educ1b = est22* educ1
gen educ2b = est23* educ2
gen educ3b = est24* educ3
gen educ4b = est25* educ4
gen educ5b = est26* educ5
gen educ6b = est27* educ6
gen educ7b = est28* educ7
gen educ8b = est29* educ8
gen educ9b = est30* educ9
gen income1b = est31* income1
gen income2b = est32* income2
gen income3b = est33* income3
gen income4b = est34* income4
gen income5b = est35* income5
gen income6b = est36* income6
gen income7b = est37* income7
gen income8b = est38* income8
gen income9b = est39* income9
gen income10b = est40* income10
gen income11b = est41* income11
gen income12b = est42* income12
gen s2b = est43*s2
gen s3b = est44*s3
gen s7b = est45*s7
gen s9b = est46*s9
gen s10b = est47*s10
gen s11b = est48*s11
gen s13b = est49*s13
gen s14b = est50*s14
gen s15b = est51*s15
gen s16b = est52*s16
gen s17b = est53*s17
gen s18b = est54*s18
gen s19b = est55*s19
gen s20b = est56*s20
gen s21b = est57*s21
gen s22b = est58*s22
gen s23b = est59*s23
gen s24b = est60*s24
gen s25b = est61*s25 
gen s26b = est62*s26
gen s27b = est63*s27
gen s28b = est64*s28
gen s29b = est65*s29
gen s30b = est66*s30
gen s31b = est67*s31
gen s32b = est68*s32
gen s35b = est69*s35
gen s36b = est70*s36
gen s37b = est71*s37
gen s38b = est72*s38
gen s39b = est73*s39
gen s40b = est74*s40
gen s41b = est75*s41
gen s42b = est76*s42
gen s43b = est77*s43
gen s44b = est78*s44
gen s46b = est79*s46
gen s47b = est80*s47
gen s48b = est81*s48
gen s49b = est82*s49
gen zmeanb = est83*zmean
gen waterb = est84*water
gen devb = est85*dev
gen forestb = est86*forest
gen shrubb = est87*shrub
gen herbb = est88*herb
gen cultb = est89*cult
gen constb = est90

preserve 
keep id est1-est90
sa "$param\nlp.dta", replace
restore

*Average Marginal Effect
generate nexpxb = exp(ptlb + msab + msa2b + shoreb + nlresb + tmeanb + tmeansqb + precipb + r1b + r2b + r3b + r4b + r5b + r8b + sexb + age1b ///
+ age2b + age3b + age4b + age5b + age6b + educ1b + educ2b + educ3b + educ4b + educ5b + educ6b + educ7b + educ8b ///
+ educ9b + income1b + income2b + income3b + income4b + income5b + income6b + income7b + income8b + income9b + income10b + income11b + income12b ///
+ s2b + s3b + s7b + s9b + s10b + s11b + s13b + s14b + s15b + s16b + s17b + s18b + s19b + s20b + s21b + s22b + s23b + s24b + s25b ///
+ s26b + s27b + s28b + s29b + s30b + s31b + s32b + s35b + s36b + s37b + s38b + s39b + s40b + s41b + s42b + s43b + s44b + s46b ///
+ s47b + s48b + s49b + zmeanb + waterb + devb + forestb + shrubb + herbb + cultb + constb)
generate nmargp = est1*nexpxb 

*Summary for predictions
sum nexpxb 
sum nmargp

*Standard Error by Delta Method
*Second derivatives:
gen n1 = est1*ptl*nexpxb + nexpxb 
gen n2 = est1* msa*nexpxb
gen n3 = est1* msa2*nexpxb
gen n4 = est1* shore*nexpxb
gen n5 = est1* nlres*nexpxb
gen n6 = est1* tmean*nexpxb
gen n7 = est1* tmeansq*nexpxb
gen n8 = est1* precip*nexpxb
gen n9 = est1* r1*nexpxb
gen n10 = est1* r2*nexpxb
gen n11 = est1* r3*nexpxb
gen n12 = est1* r4*nexpxb
gen n13 = est1* r5*nexpxb
gen n14 = est1* r8*nexpxb
gen n15 = est1* sex*nexpxb
gen n16 = est1* age1*nexpxb
gen n17 = est1* age2*nexpxb
gen n18 = est1* age3*nexpxb
gen n19 = est1* age4*nexpxb
gen n20 = est1* age5*nexpxb
gen n21 = est1* age6*nexpxb
gen n22 = est1* educ1*nexpxb
gen n23 = est1* educ2*nexpxb
gen n24 = est1* educ3*nexpxb
gen n25 = est1* educ4*nexpxb
gen n26 = est1* educ5*nexpxb
gen n27 = est1* educ6*nexpxb
gen n28 = est1* educ7*nexpxb
gen n29 = est1* educ8*nexpxb
gen n30 = est1* educ9*nexpxb
gen n31 = est1* income1*nexpxb
gen n32 = est1* income2*nexpxb
gen n33 = est1* income3*nexpxb
gen n34 = est1* income4*nexpxb
gen n35 = est1* income5*nexpxb
gen n36 = est1* income6*nexpxb
gen n37 = est1* income7*nexpxb
gen n38 = est1* income8*nexpxb
gen n39 = est1* income9*nexpxb
gen n40 = est1* income10*nexpxb
gen n41 = est1* income11*nexpxb
gen n42 = est1* income12*nexpxb
gen n43 = est1*s2*nexpxb
gen n44 = est1*s3*nexpxb
gen n45 = est1*s7*nexpxb
gen n46 = est1*s9*nexpxb
gen n47 = est1*s10*nexpxb
gen n48 = est1*s11*nexpxb
gen n49 = est1*s13*nexpxb
gen n50 = est1*s14*nexpxb
gen n51 = est1*s15*nexpxb
gen n52 = est1*s16*nexpxb
gen n53 = est1*s17*nexpxb
gen n54 = est1*s18*nexpxb
gen n55 = est1*s19*nexpxb
gen n56 = est1*s20*nexpxb
gen n57 = est1*s21*nexpxb
gen n58 = est1*s22*nexpxb
gen n59 = est1*s23*nexpxb
gen n60 = est1*s24*nexpxb
gen n61 = est1*s25*nexpxb
gen n62 = est1*s26*nexpxb
gen n63 = est1*s27*nexpxb
gen n64 = est1*s28*nexpxb
gen n65 = est1*s29*nexpxb
gen n66 = est1*s30*nexpxb
gen n67 = est1*s31*nexpxb
gen n68 = est1*s32*nexpxb
gen n69 = est1*s35*nexpxb
gen n70 = est1*s36*nexpxb
gen n71 = est1*s37*nexpxb
gen n72 = est1*s38*nexpxb
gen n73 = est1*s39*nexpxb
gen n74 = est1*s40*nexpxb
gen n75 = est1*s41*nexpxb
gen n76 = est1*s42*nexpxb
gen n77 = est1*s43*nexpxb
gen n78 = est1*s44*nexpxb
gen n79 = est1*s46*nexpxb
gen n80 = est1*s47*nexpxb
gen n81 = est1*s48*nexpxb
gen n82 = est1*s49*nexpxb
gen n83 = est1*zmean*nexpxb
gen n84 = est1*water*nexpxb
gen n85 = est1*dev*nexpxb
gen n86 = est1*forest*nexpxb
gen n87 = est1*shrub*nexpxb
gen n88 = est1*herb*nexpxb
gen n89 = est1*cult*nexpxb
gen n90 = est1*nexpxb

forvalues j=1/90{
egen nmargp`j' =mean(n`j')
}
*
forvalues k=1/90{
local nmp`k'=nmargp`k'
}
*

matrix list e(V)
matrix V = e(V)
matrix rownames V = 1-90
matrix colnames V = 1-90
matrix V1 = V[1..90,1..90]
 
matrix D_X= (`nmp1', `nmp2', `nmp3', `nmp4', `nmp5', `nmp6', `nmp7', `nmp8', `nmp9', `nmp10', `nmp11', `nmp12', `nmp13', `nmp14', `nmp15', `nmp16', `nmp17', `nmp18', `nmp19', `nmp20', ///
			 `nmp21', `nmp22', `nmp23', `nmp24', `nmp25', `nmp26', `nmp27', `nmp28', `nmp29', `nmp30', ///
			 `nmp31', `nmp32', `nmp33', `nmp34', `nmp35', `nmp36', `nmp37', `nmp38', `nmp39', `nmp40', ///
			 `nmp41', `nmp42', `nmp43', `nmp44', `nmp45', `nmp46', `nmp47', `nmp48', `nmp49', `nmp50', ///
			 `nmp51', `nmp52', `nmp53', `nmp54', `nmp55', `nmp56', `nmp57', `nmp58', `nmp59', `nmp60', ///
			 `nmp61', `nmp62', `nmp63', `nmp64', `nmp65', `nmp66', `nmp67', `nmp68', `nmp69', `nmp70', ///
			 `nmp71', `nmp72', `nmp73', `nmp74', `nmp75', `nmp76', `nmp77', `nmp78', `nmp79', `nmp80', ///
			 `nmp81', `nmp82', `nmp83', `nmp84', `nmp85', `nmp86', `nmp87', `nmp88', `nmp89', `nmp90')


matrix DVD = D_X*V1*D_X'
local nsemarg = DVD[1,1]
local nsemargse = (`nsemarg'^.5)
disp `nsemargse'

